#include "HarmonySearch.h"

#include <iostream>
#include <utility>
#include <algorithm>
using namespace std;

bool isAllZeros(const vector<bool>& vec) {
  int cnt = 0;  
  for (vector<bool>::const_iterator it = vec.begin(); it < vec.end(); ++it) {
    cnt += *it; 
  }
  return (cnt == 0);
}

double getExecTime(const SysState& sysstate, const vector<bool>& soln) {
  int n = soln.size();

  double scan_time = 0;
  for (int i = 0; i < n; ++i) {
    scan_time += soln[i] * sysstate.Ts[i];
  }

  double processing_time = 0;
  for (int i = 0; i < n; ++i)
    for (int j = i; j < n; ++j) {
      processing_time += soln[i] * soln[j] * sysstate.Tn[i][j] * sysstate.B[i][j];
  }

  return (scan_time + processing_time); 
}

bool compare_function(pair<double, int> s1, pair<double, int> s2) {
  return (s1.first < s2.first);
}

bool HarmonySearch::isInHarmonyMem(const vector<bool>& vec) {
  int hm_current_size = hmem.size();
  for (int i = 0; i < hm_current_size; ++i) {
    if (vec == hmem[i]) {
      return true;
    }
  }
  return false;
}

HarmonySearch::HarmonySearch(int hms, double hmcr, double par, double alpha, double beta, int num_files, const SysState& sysstate):
  HMS(hms), HMCR(hmcr), PAR(par), ALPHA(alpha), BETA(beta), n(num_files), sys_state(sysstate) {}

HarmonySearch::~HarmonySearch() {}

// Initializes harmony memory with HMS valid solutions.
void HarmonySearch::initHarmonyMemory() {
  for (int i = 0; i < HMS; ++i) {
    vector<bool> tmp_soln;
    while (true) {
      for (int j = 0; j < n; ++j) {
        tmp_soln.push_back(rand() % 2);
      }
      // If (1) satisfies constraint (2) is non-zero solution (3) is not currently in harmony memory, continue to generate next initial solution.
      if (!isAllZeros(tmp_soln) && !isInHarmonyMem(tmp_soln) && satisfyConstraint(tmp_soln)) {
        break;
      }
      tmp_soln.clear();
    }
    hmem.push_back(tmp_soln);
  }
  // Sorts initial solutions in terms of objective function values.
  sortHarmonyMemory();

  // Print initial best solution and associated cost.
  /*
  cout << "Initial best solution is" << endl;  
  for (int i = 0; i < n; ++i) {
    cout << hmem[0][i] << " ";
  } 
  cout << endl;
  cout << "Associated hybrid cost is" << endl;
  cout << obj_hybrid_pwt(hmem[0]) << endl;
  */
}

void HarmonySearch::sortHarmonyMemory() {
  vector<pair<double, int> > scores;
  for (int i = 0; i < HMS; ++i) {
    scores.push_back(make_pair(obj_hybrid_pwt(hmem[i]), i));
  } 
  sort(scores.begin(), scores.end(), compare_function);
  vector<vector<bool> > hmem_copy = hmem; 
  hmem.clear();
  for (int i = 0; i < HMS; ++i) {
    hmem.push_back(hmem_copy[scores[i].second]);
  }
}

bool HarmonySearch::satisfyConstraint(const vector<bool>& soln) {
  double total_scan_time = 0;
  for (int i = 0; i < n; ++i) {
    total_scan_time += soln[i] * sys_state.Ts[i];
  }
  double scan_time_bound = 0;
  for (int i = 0; i < n; ++i) {
    scan_time_bound += BETA * sys_state.Ts[i];
  }
  return (total_scan_time <= scan_time_bound);
}

double HarmonySearch::getD1(const vector<bool>& soln) {
  double exec_time = getExecTime(sys_state, soln);
  double single_file = 0;
  for (int i = 0; i < n; ++i) {
    single_file += soln[i] * sys_state.B[i][i] * (exec_time - (sys_state.Ts[i] + sys_state.Tn[i][i]));
  }
  double double_file = 0;
  for (int i = 0; i < n; ++i)
    for (int j = i + 1; j < n; ++j) {
      double_file += soln[i] * soln[j] * sys_state.B[i][j] * (exec_time - (sys_state.Ts[i] + sys_state.Ts[j] + sys_state.Tn[i][j]));
  }
  return single_file + double_file;
}

double HarmonySearch::getD2(const vector<bool>& soln) {
  double exec_time = getExecTime(sys_state, soln);
  double single_file = 0;
  for (int i = 0; i < n; ++i) {
    single_file += sys_state.L[i][i];
  }

  double double_file = 0;
  for (int i = 0; i < n; ++i)
    for (int j = i + 1; j < n; ++j) {
      double_file += sys_state.L[i][j];
  }


  return (exec_time * exec_time / 2) * (single_file + double_file);
}

double HarmonySearch::getD3(const vector<bool>& soln) {
  double exec_time = getExecTime(sys_state, soln);
  double single_file = 0;
  for (int i = 0; i < n; ++i) {
    single_file += (1 - soln[i]) * sys_state.B[i][i];
  }

  double double_file = 0;
  for (int i = 0; i < n; ++i)
    for (int j = i + 1; j < n; ++j) {
      double_file += (1 - soln[i] * soln[j]) * sys_state.B[i][j];
  }
  return exec_time * (single_file + double_file);
}

double HarmonySearch::getExpr1(const vector<bool>& soln, int i) {
  double exec_time = getExecTime(sys_state, soln);
  double nominator =
    (1 - soln[i]) * 
    ((sys_state.B[i][i] + exec_time * sys_state.L[i][i]) * (sys_state.B[i][i] + exec_time * sys_state.L[i][i] - 1) * sys_state.Tn[i][i] + exec_time * (sys_state.B[i][i] + exec_time * sys_state.L[i][i] / 2));

  double denominator = exec_time * sys_state.L[i][i] / 2 + sys_state.B[i][i];
  return nominator / denominator;
}

double HarmonySearch::getExpr2(const vector<bool>& soln, int i, int j) {
  double exec_time = getExecTime(sys_state, soln);
  double nominator =
    (1 - soln[i] * soln[j]) * 
    ((sys_state.B[i][j] + exec_time * sys_state.L[i][j]) * (sys_state.B[i][j] + exec_time * sys_state.L[i][j] - 1) * sys_state.Tn[i][j] + exec_time * (sys_state.B[i][j] + exec_time * sys_state.L[i][j] / 2));

  double denominator = exec_time * sys_state.L[i][j] / 2 + sys_state.B[i][j];
  return nominator / denominator;
}

double HarmonySearch::obj_hybrid_pwt(const vector<bool>& soln) {
  return obj_avg_pwt(soln);
//  return (obj_avg_pwt(soln) * ALPHA + obj_max_pwt(soln) * (1 - ALPHA));
}

double HarmonySearch::obj_avg_pwt(const vector<bool>& soln) {
  double D1 = getD1(soln);
  double D2 = getD2(soln);
  double D3 = getD3(soln);

  double total_num_jobs = 0;
  for (int i = 0; i < n; ++i)
    for (int j = i; j < n; ++j) {
      total_num_jobs += soln[i] * soln[j] * sys_state.B[i][j];
  }
     
  if (total_num_jobs == 0) {
    return 0;
  } else {
    return (D1 + D2 + D3) / total_num_jobs;
  }
}

// Old implementation of obj_max_pwt.
//double HarmonySearch::obj_max_pwt(const vector<bool>& soln) {
//  double max_pwt = -1.0;
//  for (int i = 0; i < n; ++i) {
//    double expr1 = getExpr1(soln, i);
//    if (expr1 > max_pwt) {
//      max_pwt = expr1;
//    }
//  }
//  for (int i = 0; i < n; ++i)
//    for (int j = i + 1; j < n; ++j) {
//      double expr2 = getExpr2(soln, i, j);
//      if (expr2 > max_pwt) {
//        max_pwt = expr2;
//      }
//  }
//  return max_pwt;
//}

// New implementation of obj_max_pwt.
double HarmonySearch::obj_max_pwt(const vector<bool>& soln) {
//  cout << "in obj_max_pwt" << endl;
  double exec_time = getExecTime(sys_state, soln);

  double max_pwt = -1.0;
  for (int i = 0; i < n; ++i) {
    for (int j = i; j < n; ++j) {
      double  tmp_max_pwt= (1 - soln[i] * soln[j]) * (sys_state.q_max[i][j]   + exec_time);  
      if (tmp_max_pwt > max_pwt) {
        max_pwt = tmp_max_pwt;
      }
    }
  } 

  double select_max_pwt = -1.0;
  for (int i = 0; i < n; ++i) {
    for (int j = i; j < n; ++j) {
			//cout << i << j << sys_state.q_max[i][j] << endl;
      double tmp_max_pwt = soln[i] * soln[j] * sys_state.q_max[i][j]; 
      if (tmp_max_pwt > select_max_pwt) {
        select_max_pwt = tmp_max_pwt;
      }
    }
  }
 
   double unselect_max_pwt = 100000000;
   for (int i = 0; i < n; ++i) {
     for (int j = i; j < n; ++j) {
       double tmp_max_pwt = (1 - soln[i] * soln[j]) * sys_state.q_max[i][j];
       if (tmp_max_pwt < unselect_max_pwt) {
         unselect_max_pwt = tmp_max_pwt;
       }
     }
   }
 
   if (select_max_pwt != 0) {
     return unselect_max_pwt / select_max_pwt;
   } else {
     //cout << "select_max_pwt is 0" << endl;
     return unselect_max_pwt;
   }
}

vector<bool> HarmonySearch::generateOneNewSoln() {
  vector<bool> new_soln(n, false) ;
  for (int i = 0; i < n; ++i) {  // for each dimension of a solution
    double hmr = rand() / RAND_MAX;
    if (hmr > HMCR) {
      // Chooses a value from harmony memory uniformly.
      new_soln[i] = hmem[rand() % HMS][i];

      // Pitch adjustment.
      double par = rand() / RAND_MAX;
      if (par > PAR) {  // Pitch adjustment with prob. PAR.
        new_soln[i] = !new_soln[i];
      }
    } else {  // Randomly chooses a value from value range {0, 1} with prob. 1 - HMCR.
      // Biased in choosing zero.
      double choose_zero = rand() / RAND_MAX;
      if (choose_zero > BETA) { 
        new_soln[i] = 0;
      } else {
        new_soln[i] = 1;
      }
    }
  }
  return new_soln;
}

vector<bool> HarmonySearch::runHarmonySearch() {
  initHarmonyMemory();

  for (int i = 0; i < NUM_ITERATIONS; ++i) {  // Termination condition: run for NUM_ITERATIONS times.
      vector<bool> new_soln = generateOneNewSoln();
      if (satisfyConstraint(new_soln) &&
          !isAllZeros(new_soln) && // avoid introducing all 0's solution into HM
          !isInHarmonyMem(new_soln)) { 
          double hybrid_score = obj_hybrid_pwt(new_soln);
          // if hybrid_score is 0, there are no jobs in the selected queues.
          if (hybrid_score != 0 && hybrid_score < obj_hybrid_pwt(hmem[HMS - 1])) {
          // If better than the worst solution in harmony memory, replaces it and then continues to next iteration.
          hmem[HMS - 1].clear();
          hmem[HMS - 1] = new_soln;
          sortHarmonyMemory();
        }
      }
  }
 
  // Creates return value from hmem[0], which is the best solution in Harmony memory.
  vector<bool> rt_vec = hmem[0];

  // Flushes harmony memory.
  hmem.clear();
  return rt_vec;
}
